#ifndef _EQUATIONSOLVER_H_
#define _EQUATIONSOLVER_H_
#include <iostream>
#include <cmath>
using namespace std;

#define infty 1e9;
const double tao = 1e-8;

//抽象函数类
class Function{
public:
    //重载运算符()来计算f()
    virtual double operator() (const double &x) const =0;

    //f()的导函数diff()
    virtual double diff(const double &x) const{
        return ((*this)(x+tao)-(*this)(x-tao))/(2*tao);};
        
};


//用来计算函数的根de父类，子类要重写求根solve()
class EquationSolver{
public:
    virtual double solve() =0;
};


//二分法求根
class BisectionSolver : public EquationSolver{
private:
    /* data */
    double a,b;
    Function &f;
    int M=10000;
    double dt=1e-7,eps=1e-8;

public:
    BisectionSolver(double a,double b,Function &f):
    f(f),a(a),b(b){};

    double solve(){
        double h,u;
        int k;
        h = b-a;
        u = f(a);
        double root,w;
        for(k=0;k<=M;k++){
            h=h/2;
            root = a+h;
            w = f(root);
            if( (h<dt && h> -dt)|| k==M) break;
            else if(w<eps && w>-eps) break;
            else if( w*u >0) a = root; 
        }
        return root;
    }
};

//牛顿迭代求根
class NewtonSolver : public EquationSolver{
private:
    /* data */
    double x0;
    Function &f;
    int M=10000;
    double eps=1e-8;

public:
    NewtonSolver(double x0,Function &f):
    f(f),x0(x0) {};

    double solve(){
        double x=x0;
        double u;
        for(int k=0;k<=M;k++){
            u = f(x);
            if(u<eps && u>-eps) break;
            x = x-u/f.diff(x);
        }
        return x;
    }
};

//割线法求根
class SecandSolver: public EquationSolver
{
private:
    /* data */
    double x0,x1;
    Function &f;
    int M=10000;
    double eps=1e-8;

public:
    SecandSolver(double x0,double x1,Function &f):
        f(f),x0(x0),x1(x1) {};

    double solve(){
        double f0= f(x0);
        double f1= f(x1);
        double s;
        for(int k=2;k<=M;k++){
            s = (x1-x0)/(f1-f0);
            x0 =x1;
            f0 =f1;
            x1 =x1 - f1*s;
            f1 = f(x1);
            if ((x1-x0)<eps && (x0-x1)<eps) break;
            if (f0<eps &&f0>-eps) break; 
        }
        return x1;
    }
};



#endif